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Abstract 

Autocatalytic reaction system with a small number of molecules is studied numerically by stochastic 
particle simulations. A novel state due to fluctuation and discreteness in molecular numbers is found, 
characterized as extinction of molecule species alternately in the autocatalytic reaction loop. Phase 
transition to this state with the change of the system size and flow is studied, while a single-molecule 
switch of the molecule distributions is reported. Relevance of the results to intracellular processes are 
briefly discussed. 

PACS numbers: 87.16.-b, 05.40.-a 



Cellular activities are supported by biochemi- 
cal reactions in a cell. To study biochemical dy- 
namic processes, rate equation for chemical reac- 
tions are often adopted for the change of chemical 
concentrations. However, the number of molecules 
in a cell is often rather small Q| , and it is not trivial 
if the rate equation approach based on the contin- 
uum limit is always justified. For example, in cell 
transduction even a single molecule can switch the 
biochemical state of a cell In our visual sys- 
tem, a single photon in retina is amplified to a 
macroscopic level HJ . 

Of course, fluctuations due to a finite number 
of molecules is discussed by stochastic differen- 
tial equation (SDE) adding a noise term to the 
rate equation for the concentration |^. This 
noise term sometimes introduces a non-trivial ef- 
fect, as discussed as noise-induced phase transition 
|3|, noise- induced order 0, stochastic resonance 
|B|, and so forth. Still, these studies assume that 
the average dynamics are governed by the contin- 
uum limit, and the noise term is added as a per- 
turbation to it. 

In a cell, often the number of some molecules 
is very small, and may go down very close to or 
equal to 0. In this case, the change of the number 
between zero and nonzero, together with the fluc- 
tuations may cause a drastic effect that cannot be 
treated by SDE. Possibility of some order differ- 
ent from macroscopic dissipative structure is also 
discussed by Mikhailov and Hess |l^ (see also 



Ref . |llj] ) . Here we present a simple example with 
a phenomenon intrinsic to a system with a small 
number of molecules where both the fluctuations 
and digitality('0/l') are essential. 

In nonlinear dynamics, drastic effect of a sin- 
gle molecule may be expected if a small change is 
amplified. Indeed, autocatalytic reaction widely 
seen in a cell, provides a candidate for such am- 
plification jl^, Here we consider the simplest 
example of autocatalytic reaction networks (loops) 
with a non-trivial finite- number effect. With a cell 
in mind, we consider reaction of molecules in a 
container, contacted with a reservoir of molecules. 
The autocatalytic reaction loop is Xi -I- ATj+i 
2Xi+i;i = 1, • • • , fc; Xk+i = Xi 
within a container. Through the contact with a 
reservoir, each molecule Xi diffuses in and out. 

Assuming that the chemicals are well stirred in 
the container, our system is characterized by the 
number of molecules Ni of the chemical Xi in the 
container with the volume V In the contin- 

uum limit with a large number of molecules, the 
evolution of concentrations Xi = Ni/V is repre- 
sented by 

dxi/dt = riXj_iXj -Ti+iXjXj+i + A(si - Xi) (1) 

where is the reaction rate, Di is the diffusion 
rate across the surface of the container, and Si is 
the concentration of the molecule in the reservoir. 

For simplicity, we consider the case = r, 
Di — D, and s,; = s for all i, while the phenom- 
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ena to be presented here will persist by dropping 
this condition. With this homogeneous parameter 
case, the above equation has a unique attractor, a 
stable fixed point solution with Xi = s. The Jacobi 
matrix around this fixed point solution has a com- 
plex eigenvalue, and the fluctuations around the 
fixed point relax with the frequency tOp = rs/ir. 
In the present paper we mainly discuss the case 
with fc = 4, since it is the minimal number to see 
the new phase to be presented. 

If the number of molecules is finite but large, 
the reaction dynamics can be replaced by Langevin 
equation by adding a noise term to eq. ([|). In 
this case, the concentration Xi fluctuates around 
the fixed point, with the dynamics of a component 
of the frequency cOp. No remarkable change is ob- 
served with the increase of the noise strength, that 
corresponds to the decrease of the total number of 
molecules. 

To study if there is a phenomenon that is out- 
side of this SDE approach, we have directly sim- 
ulated the above autocatalytic reaction model, by 
colliding molecules stochastically. Taking randomly 
a pair of particles and examining if they can react 
or not, we have made the reaction with the prob- 
ability proportional to r. On the other hand, the 
diffusion out to the reservoir is taken account of 
by randomly sampling molecules and probabilis- 
tically removing them with in proportion to the 
diffusion rate D, while the fiow to the container 
is also carried out stochastically in proportion to 
s, D and V Technically, we divide time into 
time interval St for computation, where one pair 
for the reaction, and single molecules for diffusion 
in and out are checked. The state of the container 
is always updated when a reaction or a flow of a 
molecule has occurred. The reaction Xi -l-X^+i 
2Xi+i is made with the probability Piii{t,t + St) = 
rx,{t)x^+i(t)Vdt = rN,(t)N,+i{t)V-^dt within the 
step 6t. A molecule diffuses out with the proba- 
bility Poi = DVxi{t) = DNi{t), and flows in with 
Pji = DVs. We choose 6t small enough so that the 
numerical result is insensitive with the further de- 
crease of St. By decreasing Vs, we can control the 
average number of molecules in the container, and 
discuss the effect of a finite number of molecules, 
since the average of the total number of molecules 
Ntot is around the order of AVs jl^ . On the other 
hand, the 'discreteness' in the diffusion is clearer 
as the diffusion rate D is decreased. We set r = 1 
and s = 1, without loss of generality {rs/D and 
sV are the only relevant parameters of the model 
by properly scaling the time t). 

First, our numerical results agree with those 
obtained by the corresponding Langevin equation 
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Figure 1: (color) Time series of the number of 
molecules Ni{t), for D = 1/256, V = 32. Either 1- 
3 or 2-4 rich state is stabilized. Successive switches 
appear between A^i > 7V3 and N3 > Ni states with 
N2, Ni « 0. Here a switch from 1-3 rich to 2-4 rich 
state occurs around t ~ 5500. 

if D and V are not too small. As the volume V 
(and accordingly Ntot) is decreased, however, we 
have found a new state whose correspondent does 
not exist in the continuum limit. An example of 
the time series is plotted in Fig. |l|, where we note 
a novel state with Ni, ^ 1 and N2, N4 oi or 
A^2, A'4 > 1 and Ni,N3 w 0. To characterize this 
state quantitatively, we have measured the prob- 
ability distribution of z = zi -f 2:3 — {x2 + 2:4). 
Since the solution of the continuum limit is Xi ~ 
s{= 1) for all i, this distribution has a sharp peak 
around 0, with a Gaussian form approximately, 
when Ntot'is large enough. As shown in Fig. the 
distribution starts to have double peaks around 
±4, as V is decreased. With the decrease of V (i.e., 
Ntot), these double peaks first sharpen, and then 
get broader with the further decrease due to too 
large fluctuation of a system with a small number 
of molecules. Hence the new state with switches 
between 1-3 rich and 2-4 rich temporal domains 
is a characteristic phenomenon that appears only 
within some range of a small number of molecules 

The stability of this state is understood as fol- 
lows. Consider the case with 1-3 rich and N2 = 
N4 = 0. When one (or few) X2 molecules flow 
in, N2 increases, due to the autocatalytic reaction. 
Then X3 is amplified, and since N2 is not large, 
N2 soon comes back to again. In short, switch 
from (7Vi,0, A^3,0) to (iVi - A,0,iV3 + A + 1,0) 
occurs with some A, but the 1-3 rich state itself 
is maintained. In the same manner, this state is 
stable against the flow of X4. The 1-3 rich state is 
maintained unless either iVi or A'3 is close or equal 
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D = 1/64 
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Figure 2: (color) The probability distribution of 
z = (x\ + x-i) — (xi + X4), sampled over 2.1 - 
5.2 X 10*^ steps. D = 1/64. For V > 128, z has a 
distribution around 0, corresponding to the fixed 
point state Xi = s{~ 1). For V < 32, the dis- 
tribution has double peaks around z w ±4, cor- 
responding to the state iVi,iV3 > 7V2,^4(~ 0) 
or the other way round. The double-peak distri- 
bution is sharpest around V = 16, and with the 
further decrease of V, the distribution is broader 
due to finite-size fluctuations. 

to 0, and both X2 and X4 molecules flow in within 
the switch time. Hence the 1-3 rich state (as well 
as 2-4 rich state, of course) is stable as long as the 
flow rate is small enough. 

Within a temporal domain of 1-3 rich state, 
switches occur to change from (A^i, N3) ^ {N[, N^) 
In Fig. ^, we have plotted the probability density 
for the switch from A^i N[ when a single X2 
molecule flows in, amplified, and N2 comes back to 
0, by fixing Ni + N3 = Nini at 256 initially (We 
assume no more flow. Hence -|- = N^i -1-1). 
The peak around N{ ^ Ni + 1 means the reaction 
from N2 to before the amplification, while an- 
other peak around N[ ^ = Nini — Ni shows 
the conversion of the numbers through the ampli- 
fication of X2 molecules. Indeed, each temporal 
domain of the 1-3 rich state consists of successive 
switches of (NijN^) (iV3,iVi), as shown in 

Fig. 1^. Since molecules diffuse out or in randomly 
besides this switch, the difference between A^i and 
N3 is tended to decrease. On the other hand, each 
1-3 rich state, when formed, has imbalance be- 
tween A^i and iVa, i.e., iVi > N3 or iVi < N3, 
since, as in Fig. the state is attracted from al- 
ternate amplification of Xi, where only one type i 
of molecules has Ni ^ 1 and for others. How- 
ever, the destruction of the 1-3 rich state is easier 
if Ni :s> N3 or iVi ^ iVs, as mentioned. Roughly 
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Figure 3: Probability density for the switch from 
{Ni^Ns) to {N[,N^) when a single X2 molecule is 
injected into the system. Ni + N3 = Nini is fixed 
at 256 initially. There is no more fiow and N4 is 
always kept at 0, so that the switch is completed 
when N2 comes back to 0, and A''^ -I- A^g = Nini + 1- 
Probability to take N[ is plotted against initial iVi . 

speaking, each 1-3 rich state starts with a large 
imbalance between Ni and N3 , and continues over 
a long time span, if the switch and diffusion lead to 
iVi « iVg, and is destroyed when the large imbal- 
ance is restored. Indeed, we have plotted the dis- 
tribution of y = xi—X3 + X2 — X4^, to see the imbal- 
ance for each 1-3 rich or 2-4 rich domain. This dis- 
tribution shows double peaks clearly around y « 
±2.8,i.e.,{Ni,N3) « (3.4^, 0.6^), (0.6^, 3.4^). 

Let us now discuss the condition to have the 
1-3 or 2-4 rich state. First, the total number of 
molecules should be small enough so that the fiuc- 
tuation from the state Ni « Nj (for Vi,j) may 
reach the state with Ni « 0. On the other hand, 
if the total number is too small, even iVi or N3 for 
the 1-3 rich state may approach easily, and the 
state is easily destabilized. Hence the alternately 
rich state is stabilized only within some range of 
V. 

Note also that our system has conserved quan- 
tities Ni (and logxi in the continuum limit), 
if D is set at 0. Hence, as the diffusion rate gets 
smaller, some characteristics of the initial popu- 
lation are maintained over long time. Once the 
above 1-3 (or 2-4) rich state is formed, it is more 
difficult to be destabilized if D is small. In Fig. 
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Figure 4: The rate of the residence at 1-3 (or 2-4) 
rich state over the whole temporal domain, plotted 
against DV jlSj. Here, the residence rate is com- 
puted as follows. As long as -/V2 > and N4 > 
are not satisfied simultaneously, over a given time 
interval (8.0; 2.5 times as long as the period of 
the oscillation around the fixed point at contin- 
uum limit), it is counted as the 1-3 rich state (2-4 
rich state is defined in the same way) . The res- 
idence rate is computed as the ratio of the fraction 
of the time intervals of 1-3 or 2-4 rich state to the 
whole interval. For small D, this switching state 
is observed even for a large number of molecules, 
say = 4 X lO'^, for V ^ 10^ at D ^ 10"^. 



we have plotted the rate of the residence at 1-3 
(or 2-4) rich state over the whole temporal domain, 
with the change of V. Roughly speaking, the state 
appears for DV < 1 [Q, while for too small V 
(e.g., F < 4), it is again destabilized by fluctua- 
tions. Although the range of the 1-3 rich state is 
larger for small D, the necessary time to approach 
it increases linearly with V. Hence it would be fair 
to state that properly small number of molecules 
is necessary to have the present state. 

To sum up, we have discovered a novel state 
in reaction dynamics intrinsic to a small number 
of molecules. This state is characterized by alter- 
nately vanishing chemicals within an autocatalytic 
loop, and switches by a flow of single molecules 
^ . Hence, this state generally appears for a sys- 
tem with an autocatalytic loop consisting of any 
even number of elements. With the increase of k, 
however, the globally alternating state all over the 
loop is more difficult to be reached. In this case, 
locally alternating states are often formed with the 
decrease of the system size (e.g., '2-4-6-8 rich' and 
'11-13-15 rich' states for k = 16). This local order 
is more vulnerable to the fiow of molecules than 
the global order for the k — A loop. 

On the other hand, for fc = 3, two of the chem- 



ical species start to vanish for small V^ since any 
pair of different chemical species can react so that 
one chemical species is quickly absorbed into the 
other. This state of single chemical species, how- 
ever, is not stable by a fiow of a single molecule. 
Indeed, no clear 'phase transition' is observed with 
the decrease of V. 

Although in the present Letter we have studied 
the case with Sj = s, we have also confirmed that 
the present state with alternately vanishing chem- 
ical species is generally stabilized for small V, even 
if Si or Vi or Di are not identical. 

Last, we make a remark about the signal trans- 
duction in a cell. In a cell, often the number of 
molecules is small, and the cellular states often 
switch by a stimulus of a single molecule ^ . Fur- 
thermore, signal transduction pathways generally 
include autocatalytic reactions. In this sense, the 
present stabilization of the alternately rich state as 
well as a single-molecule switch may be relevant 
to cellular dynamics. Of course, one may won- 
der that the present mechanism is too 'stochastic'. 
Then, use of both the present mechanism and ro- 
bustness by dynamical systems ||2C , 21 may be im- 
portant. Indeed, we have made some preliminary 
simulations of complex reaction networks. Often, 
we have found the transition to a new state at 
a small number of molecules, when the network 
includes the autocatalytic loop of 4 chemicals as 
studied here Hence the state presented here 
is not restricted to this specific reaction network, 
but is observed in a class of autocatalytic reaction 
network. Furthermore switches between different 
dynamic states (limit cycles or chaos) are possible 
when the number of some molecules (that are not 
directly responsible to the switch) is large enough. 
The 'switch of dynamical systems' by the present 
few-number-molecule mechanism will be an impor- 
tant topic to be pursued in future. 
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